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ABSTRACT 


—  This  paper  gives  a  general  method  for  the  stable  evaluation  of 
multivariate  B- splines.  The  problem  of  evaluation  along  mesh  boundaries  is 
discussed  in  detail.  Several  examples  are  presented  to  demonstrate  the 


effectiveness  of  the  method  for  arbitrary  B-splines 
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SIGNIFICANCE  AND  EXPLANATION 


■  y  For  on*  variable,  the  problem  of  stably  evaluating  B-splines  via  their 
recurrence  relations  is  well  understood.  For  multivariate  B-splines,  however 


the  gecsMtry  becomes  more  complex,  and  it  becomes  quite  difficult  to  implemen 


the  recurrence  relations  in  such  a  way  that  one  ends  up  with  a  robust 
evaluation  method.  .  This  paper  presents  a  method  which  guarantees  the  stable 


evaluation  of  all  smooth  multivariate  simplex  B-splines. 
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Nrtlwnow,  if  *  =  with  =  0,  thin 


Proafe  The  proof  of  thio  theorem  hot  bow  |hio  in  various  form*  bp  »  number  of  pooplo, 
mdudmg  M iccbeffi  (1979),  BoUig  [1980],  Rakopua  [1980],  and  de  Boor  end  Hoffig  [1982]. 

Grow  Theorem  1,  it  should  now  bo  possible  to  evaluate  a  spline,  b  thi^  theorem,  all  the 


iaductlvebr  write  agy  opine  aa  a 


This  is,  Ibr  a  specific  collection  of  pointa  a  straightforward  teak,  the  difficult  arieee 

when  one  wishes  to  write  a  coanpnter  code  to  evelnite  B  sphnae  taring  (B)  for  erfctrsrycoUectione 
of  potato.  The  problem,  of  comae,  is  to  feed  the  fe  so  that  s  =  and  *  1  are 

sathfeed  hi  the  general  ewe,  fat  each  a  wap  the*  neither  accuracy  nor  efficiency  sniff.  Accuracy 
becomes  an  feme  whenever  one  is  tawed  to  nes  *  of  opposite  Mg*.  Since  £**  a,  *  1,  this  maws 
that  one  aboald  imist  oh  having  o< >  0,  far  all  •  *  0,  ...n.  Efficiency  dictates  that  as  many  of  the 
as  possible  be  aero.  Thm,  one  wishes  to  feed  a  eolation  to  the  following  problem: 


jn 

x>«<=* 


w 

2>=1 

<*o 

a<5  0  ,  1  =  0,...,** 

Tbs  difficulty  with  this  Bee  in  the  non  negativity  constraints.  Efficiency  demands  that  at  most 
mil  of  the  a<  be  non-soro,  yet  arbitrarily  setting  tbs  remaining  n-m+1  of  the  a<  to  aero  provides 
no  guaranty  that  tbs  noo-negativity  constraints  will  be  satisfied. 

Example  3:  Taking  again  the  bivariate  B*spUne  giveo  by  the  points  to,  given  in  Example 

I,  and  the  point  (*,p)  =  ( j,  |),  one  might  wish  to  compute  A#((c,g)|t0,...,ts).  Setting  a*  =  a4  = 

as  *  0,  one  gets 

3 

a°*  4 

1 

0|  *  ~4 
1 

«*  *  2* 


2 


which  is  not  s  solution  to  (5). 


Fortunately,  then  is  a  well-established  procedure  for  handling  problems  of  a  similar  natum. 
Except  for  the  absense  of  an  objective  function,  (5)  is  of  the  same  form  as  a  linear  programming 
problem.  Such  problems  have  been  studied  in  gnat  detail  over  the  yean,  and  many  ways  of  com¬ 
puting  their  solutions  an  known.  Thus,  an  approach  which  one  might  consider  is  the  introduction 
of  a  "dummy”  objective  function  to  convert  (5)  to  a  linear  programming  problem  which  may  then 
be  solved  by  standard  techniques. 

In  practice,  this  seems  to  work  quite  well,  in  the  numerical  experiments  performed  to  date,  the 
simplex  method,  implemented  with  the  help  of  the  Tucker  tableau  (described  in  the  Appendix), 
has  performed  admirably.  One  introduces  the  objective  function  0  (for  masons  to  be  made  clear 
shortly)  and  tmats  (5)  as  a  linear  programming  problem;  i.e.,  one  considen  the  following  problem, 
equivalent  to  (5): 

max  0 


subject  to 


= * 

i=0 


1 


(«) 


Qi  >  0  ,  »  =  0,...,n 


Since  the  objective  function  is  constant,  solving  this  problem  amounts  to  finding  a  feasible  point. 
In  other  words,  it  amounts  to  finding  a  solution  to  (S). 

In  solving  (6)  via  the  method  outlined  in  the  preceding  paragraph,  one  can  see  why  the 
particular  choice  of  an  objective  function  is,  in  some  sense,  optimal.  In  the  Tucker  tableau,  the 
pivot  rules  an  such  that  a  row  of  all  sens  will  never  change.  Thus,  the  physical  storage  of  the 
row  corresponding  to  the  objective  function  in  memory  is  unnecessary.  Furthermore,  the  dual  of 
the  problem  is  given  by 


min  (*,«)  +  » 

subject  to  (/**<,«) +  e>0  ,  »  =  0, ...,«. 


This  problem  is  dual  feasible  in  that  u  =  0,  v  =  0  satisfies  the  constraints  for  this  problem. 
This  situation  arises  because  the  objective  function  of  the  primal  problem  (6)  has  non-positive 
coefficients.  This  means  that  (6)  is  “easy"  to  solve  because  one  can  dispense  with  phase  1  of  the 
simplex  method  and  use  instead  the  dual  simplex  method. 

Example  4:  Consider  again  Example  2.  Setting  up  this  problem  as  a  linear  programming 
problem  leads  to  the  tableau 


Oo 

on 

or* 

01$ 

«4 

1 

r0  l  -1 

-1 

0 

1 

1 

0 

_1 

s 

n  0 

-1 

-1 

0 

1 

1 

_1 

4 

ra  l-l 

-1 

-1 

-1 

-I 

-1 

-1 

where  the  variables  ro,  rj,  and  ra  are  so-called  "slack”  variables.  Since  (6)  is  made  up  of  equality 
constraints,  one  first  pivots  r0,  rlf  and  rs  to  the  top  of  the  tableau  and,  once  this  is  done,  deletes 
the  columns  corresponding  to  them.  After  all,  the  variables  along  the  top  are  assumed  to  be 


zero,  and  the  deletion  of  a  column  corresponding  to  such  a  variable  merely  makes  this  condition 
permanent.  Thus,  one  first  exchanges  r0  and  Qo  to  get 

ro  oj  O]  as  o«  o$  1 


o0/ -1  1  0-1-1  0  \ 

r,  I  0  -1  -1  0  1  1 

rt  \-l  0  -1  -2  -2  -1 


one  deletes  the  first  column  and  exchanges  n  and  a  1  to  get 
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Of 

Ol 

04 

Of 

00  /  1 

-1 

-1 

0 

1 

Ol  -1 

1 

0 

-1 

-1 

r,\0 

-1 

-2 

-2 

-I 

Lastly,  one  deletes  the  column  corresponding  to  rj,  exchanges  rj  and  as,  and  deletes  the  resulting 
column  corresponding  to  rs  to  obtain 


Os 

o4 

as 

QO  *  1 

2 

2 

aj  -2 

-3 

-2 

“a  2 

2 

1 

This  is  not  (dual)  optimal,  so  one  must  exchange  as  with  some  column  whose  entry  in  Os’s  row  is 
negative.  Since  ill  columns  satisfy  this,  the  first  is  chosen,  and  oi  and  a%  are  exchanged  to  reveal 

Oi  04  0(  1 

“0  (  \  \  in 

M  I  *  i 

o,  V  1  -1  -1  \) 

This  is  a  solution  to  (6),  and  therefore  to  (5).  Note  that  although  not  explicitly  required,  only  3, 
or  m  +  1  of  the  o »  are  non-zero.  This  occurs  because  (6)  only  had  3  constraints  (other  than  the 
non-negativity  constraints),  and  solutions  of  a  linear  programming  problem  must  satisfy  a  com¬ 
plementarity  condition;  that  is,  the  only  variables  which  can  be  non- zero  are  those  corresponding 
to  tight  constraints.  Since  r0,  rj,  and  rj  were  forced  to  be  zero,  there  can  be  at  most  3  non-zero 
a,.  Thus,  the  linear  programming  approach  implicitly  takes  care  of  the  efficiency  issue  discussed 
above. 

Once  one  has  solved  (6)  by  this  approach,  one  can  easily  evaluate  the  spline  using  (4),  assuming 
that  the  values  of  the  lower  order  splines  which  occur  on  the  right  hand  side  of  the  equality  are 
known.  In  general,  this  is  not  the  case,  but  one  can  reapply  the  technique  to  each  of  the  splines 
appearing  on  the  right  hand  side  of  (4).  This  process  may  be  carried  out  inductively  until  one  can 
finally  express  the  value  of  the  desired  B-spline  at  the  desired  point  in  terms  of  piecewise  constant 
functions  at  that  point. 

The  l\icker  tableau  makes  this  inductive  process  extremely  efficient.  The  tableau  which  solves 
(6)  may  be  used  to  solve  the  resulting  subproblems.  One  can  view  each  of  the  subproblems  as 
being  just  (5)  with  the  additional  constraint  a,  =  0.  Thus,  to  solve  a  subproblem,  one  can  take 
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the  tableau  which  solves  the  main  problem,  pivot  a,  to  the  top  of  the  tableau,  delete  the  column 
corresponding  to  it,  and  then  optimize.  This  will  yield  a  solution  to  the  subproblem  in  short  order, 
often  only  one  pivot.  These  solution  tableaus  for  the  subproblems  may  then  be  used  to  obtain 
cheap  solutions  to  the  sub-subproblems,  etc.  This  exploitation  of  similarity  among  the  various 
linear  programming  problems  which  one  must  solve  saves  an  enormous  amount  of  work. 

At  first  glance,  it  would  seem  that  the  difficulties  in  computing  the  value  of  a  multivariate  B- 
spline  have  been  overcome.  Unfortunately,  one  of  the  more  persistent  of  the  problems  has  yet  to  be 
overcome.  The  piecewise  constant  functions  which  one  ultimately  ends  up  with  have  discontinuities 
along  certain  boundaries,  namely  along  the  grid  lines.  A  grid  line  is  a  set  consisting  of  the  convex 
hull  of  m  or  fewer  points  taken  from  the  set  {Pto,  •••,  Ptn}-  A  point  x  is  said  to  lie  on  a  grid  line 
if  it  is  a  member  of  some  such  set.  Whenever  one  wishes  to  evaluate  a  spline  at  a  point  lying  on 
such  a  grid  line,  one  runs  the  risk  of  computing  it  improperly. 

Example  3:  Suppose  one  wishes  to  evaluate  A/(0,0|to, ti,*a,ts),  where  to  —  (M>0)»  = 

(-1,1,0),  t ,  =  (-1,-1,0),  and  ta  =  (1,  — 1,1).  Then  (0,0)  =  \Pt0  +  \Ptit  and  therefore 

M(0,0|t0,ti,tj,ts)  =  ^M(0,0|t1,ta,ts)+  ^M(0,0|to,ti,ts). 

But  M(*|ti,ta,ts)  and  M(«|t0,ti,ts)  are  discontinuous  at  (0,0),  so  it  is  unclear  whether  to  choose 
the  interior  or  exterior  limits  as  values  for  these  splines.  If  interior  limits  are  chosen,  the  computed 
value  of  the  spline  will  be  twice  as  great  as  the  actual  value.  If  exterior  limits  are  chosen,  the  value 
will  be  zero,  and  this  is  obviously  also  incorrect. 

There  are  many  ways  of  attempting  to  circumvent  this  problem,  but  nearly  all  fail  in  some 
way,  especially  when  one  takes  into  account  the  inexact  nature  of  the  arithmetic  performed  by  the 
computer.  In  general,  there  seems  to  be  no  reasonable  way  to  handle  this  problem,  but  for  smooth 
splines,  there  are  a  few  things  one  might  try. 

The  obvious  approach  is  to  prohibit  one  from  evaluating  a  spline  on  the  grid  lines.  This  is 
certainly  the  most  sure-fire  answer,  and  it  is  also  a  simple  enough  scheme  to  be  easily  implemented. 
All  one  need  do  whenever  one  finds  that  he  is  on  a  grid  line  is  to  move  the  point  in  some  direction 
by  e  and  try  again.  Higher  order  splines  are,  in  general,  continuous,  so  this  small  change  in  the 
location  of  the  point  will  make  a  very  small  change  in  the  value  of  the  spline.  Unfortunately,  this 
must  be  done  by  hand,  since  the  computer  is  unable  to  tell  when  a  point  is  actually  “mi”  a  grid 
line;  the  best  it  can  usually  do  is  to  tell  when  a  point  is  “near”  a  grid  line.  As  the  number  of 
variables  increases,  the  structure  of  the  grid  lines  becomes  increasingly  complex.  As  a  result,  it 
becomes  increasingly  difficult  to  avoid  computing  the  value  of  the  spline  there.  For  example,  in 
one  variable,  the  grid  lines  consist  only  of  the  knots,  while  in  two  variables,  the  grid  lines  consist 
of  the  knots  as  well  as  the  line  segments  joining  the  knots  (see  Figure  1). 

In  one  variable,  the  problem  isn’t  so  terrible.  Typically,  one  decides  that  all  the  piecewise 
constant  splines  are  either  continuous  from  the  left  or  continuous  from  the  right.  Then,  when  one 
needs  to  evaluate  at  a  knot,  one  sets  the  value  of  the  spline  to  zero  if  it  is  the  left  knot,  for  example, 
and  non-zero  if  it  is  the  right  knot.  Thus,  for  one  variable  at  least,  little  needs  to  be  done  to  clear 
up  this  nuisance. 

Conveniently,  such  a  plan  of  attack  generalizes  to  more  than  one  variable.  One  merely  chooses 
(somewhat  arbitrarily)  some  direction  in  R"*  and  evaluates  piecewise  constant  splines  according 
to  the  following  rule:  If  *  is  in  the  interior  of  the  region  of  support,  or  if  x  is  on  the  boundary 
of  the  region  of  support  and  the  arbitrarily  chosen  direction  points  into  the  interior,  the  value 
of  the  piecewise  constant  spline  shall  be  used;  in  all  other  cases  the  value  of  the  spline  shall  be 
0.  In  theory,  this  rule  eliminates  the  difficulty.  Unfortunately,  because  of  roundoff  error,  one 
can  have  both  situations  occurring,  and  the  ambiguity  about  what  to  do  persists.  Thus,  for 


multivariate  splines  an  alternate  approach  must  be  taken.  Its  successful  implementation  depends 
on  the  following  theorem. 

Theorem  2:  Let  to,—  »*»-n  be  a  collection  of  points  in  general  position  in  R".  Let  A  := 
l*o,—  ,*»+:],  end  let  Ai  :=  (to,— Then  for  all  *  €  A,  with  x  not  on  any  of  the 
grid  lines,  x  lies  in  exactly  two  of  the  it,. 

Proof:  The  statement  x  6  A  is  equivalent  to  the  statement  that  there  exists  a  solution  to  the 
problem: 


»+i 


»= o 

n+I 


* 

«=0 

a<  >  0,  »  =  0,  ...n+  1. 

Ignoring  for  the  moment  the  inequality  constraint,  one  can  view  (7)  as  a  linear  system,  namely 

Ta  =  /?, 


(7) 


where  Tisann  +  lxn  +  2  matrix,  a  is  the  vector  whose  individual  components  are  the  a^,  and  p 
is  the  vector  obtained  by  adding  the  component  1  to  the  end  of  x.  T  clearly  has  rank  n  + 1,  since 
the  (,  are  in  general  position.  Thus,  this  linear  system  has  a  one-parameter  family  of  solutions,  say 
e(i)  :=  y  +  w,  where  y,t€  Rn+1  and  «  €  R.  Now  one  can  consider  the  inequality  constraints, 
<*•  >  0, »  =  0,  ...,n  + 1.  This  is  equivalent  to  yt  +  sz,  >  0,  *  =  0,...,n  + 1.  Taken  together,  all  these 
conditions  define  some  interval  5  :=  [»_,«+]  in  which  s  must  lie  in  order  for  a(s)  to  satisfy  the 
inequality  constraints.  Since  x  €  A,  it  is  clear  that  S  is  non-empty.  Furthermore,  it  is  clear  that 
fj  =  l  and  Y*=o  *»  =  0,  for  a •  ~  * » independent  of  s.  Since  the  solution  cannot  be 
unique,  at  least  one  of  the  *,  is  non-sero.  But  <«  —  0,  so  there  must  be  at  least  two  of  the  z, 
non-sero  and  of  opposite  sign.  When  Zi  is  positive,  one  gets  a  lower  bound  for  s,  while  z<  negative 
gives  an  upper  bound  for  z.  Hence,  S  is  a  finite  interval.  Since  a(s)  is  a  continuous  function  of 
s,  no  components  of  a  can  have  sign  changes  in  S.  Furthermore,  a  outside  of  S  means  that  one 
or  more  components  of  a(«)  are  negative  there,  hence  must  change  sign  on  the  boundary  of  S. 
Suppose  <**(#_)  =  0  and  Oj(s+)  =  0.  Since  x  does  not  lie  on  a  grid  line,  <*,(*-)  is  the  only  zero 
component  of  a  at  s_ .  Similarly,  ay(«+)  is  the  only  zero  component  of  a  at  But  this  says  that 
x  €  Ai  and  x  €  Ar  Since  a(«)  is  affine,  no  other  solutions  with  a  zero  component  are  possible. 
Thus,  x  lies  in  exactly  two  of  the  Ai.  This  proves  the  theorem. 

With  this  theorem,  one  can  now  correctly  evaluate  continuous  piecewise  linear  B-splines. 
Suppose  one  wishes  to  evaluate  M(x|to,— where  Pt0,...,  Pt„+\  are  points  in  general  position 
in  Rn.  Then,  after  one  solves  (5)  to  get 

*»+l 

A/(x|to,  —  ,tn+i)  =  (n  +  1)  )  ]  a,M(x\t0,  —  >*n+i)> 

•=0 

it  is  dear  that  if  x  does  not  lie  on  a  grid  line,  then  all  but  one  of  the  Af(xjfo,  ...,t  — ,t„+i) 

will  be  zero.  This  is  an  immediate  consequence  of  Theorem  2.  If  x  is  on  a  grid  line,  however, 
one  can  still  impose  this  condition.  This  trick  forces  evaluation  on  the  grid  lines  to  behave  just 
like  evaluation  off  the  grid  lines.  One  must  be  careful,  however,  to  choose  the  correct  piecewise 
constant  spline  to  be  non- zero. 


ft  CoMi dar  again  Example  4.  One  impoeea  the  condition  that  exactly  one  of 
and  kf(e|to,ts,t*)  can  be  non-sero  at  (0,0).  la  tide  caae,  it  doesn't  matter  what  one 
>  might  aHempt  to  solve  this  problem  numerically  with  the  following  reeult: 


A*(0,0Mi»*>»*s)=  |M(0,0jtI,t,,t,)  +  10'«Af(0,0|t0,t,,t,)+  ?A#(0,0|to,t„t,). 

One  dearly  cannot  choose  M(e|t0,  *>.*>)  as  the  only  non-iero  spline  and  expect  to  get  reasonable 


This  indicates  tfca&one  must  be  careful  in  the  selection  of  the  non-sero  spline.  A  good  method 
is  to  choose  from  all  possible  splines  the  one  which  has  the  largest  coefficient.  This  will  elimhiste 
the  kind  of  numeric  nuisance  which  occurs  in  Example  5. 

However,  it  must  be  pointed  out  diet  this  approach  only  works  if  the  linear  spline  is  con¬ 
tinuous.  If  it  is  discontinuous,  one  can  evaluate  the  spline  stably  everywhere  except  along  the 
discontinuity,  where  numeric  noise  makes  the  exact  location  of  the  discontinuity  impossible  to  cal¬ 
culate.  fortunately,  lor  smooth  splines  at  least,  this  never  happens,  and  one  needn’t  be  concerned 
with  it. 

Given  that  one  can  evaluate  continuous  linear  splines  stably  everywhere,  one  can  evaluate 
smooth  higher  degree  splines  stably  everywhere.  Instead  of  expressing  the  higher  degree  spline  as 
a  linear  combination  of  constant  splines,  one  stops  one  level  sooner  and  expresses  it  as  a  linear 
combination  of  linear  splines,  each  of  which  is  continuous  and  can  be  evaluated  stably  by  the 
method  just  ^ascribed.  One  no  longer  worries  about  extraneous  terms  resulting  from  grid  line 
effects  because  the  linear  splines  are  continuous. 

Using  this  technique,  many  different  splines  have  bam  computed  in  the  bivariate  case.  Some 
of  these  appear  on  the  foUdwring  pages,  ia  order  to  produce  these  graphs,  the  mesh  was  deliberately 
chosen  so  that  many  evaluations  along  grid  lines  were  necessary.  The  knots  used  for  each  spline 
an  given  at  the  bottom  of  each  page. 


Quadratic  Bivariate  B-spline 


$ 

* 

& 


Cubic  Bivariate  B-spline 


(1,0), (1,1), (0,1),(-1,0),(- 1,-1), (0,-1) 
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Appendix:  The  Tucker  Tablean 


Consider  the  folkwing  linear  programming  problem: 

min  x  :=  (e,x)  -  d 

subject  to  r  :=  b  -  Ax  >  0  (Al) 

*<>0,  *  = 

Each  constraint,  r*  >  0  or  x«  >  0,  defines  a  half-space  in  which  the  solution  must  lie.  Since  one 
must  satisfy  all  the  constraints,  the  solution  must  lie  in  the  intersection  of  these  various  h airspaces, 
i.e.,  x  must  lie  in  some  polyhedral  region  in  R",  the  feasible  region.  The  functional  which  one 
wishes  to  minimise  over  this  region  is  linear,  so  the  solution  must  lie  at  an  extreme  point  of  the 
feasible  region,  sometimes  called  a  vertex.  A  vertex  is,  in  general,  determined  by  specifying  n  of 
the  bounding  hyperplanes  to  which  it  belongs.  This  means  that  a  vertex  is  determined  by  setting 
n  of  the  numbers  rj,  r*, ...,  rmt  x\,  **, ...,  zn  to  sero.  These  are  the  non-baafe  variables  for  that 
vertex;  call  them  £(,  (j,  ....  The  remaining  variables  are  the  basic  variables  for  that  vertex, 
and  they  will  be  denoted  by  Pi,  ps, ...,  pm.  They  are  related  by  the  equation 


P  =  0-Tt,  (A2) 

which  is  obtained  from  the  equation 

r  =  b-  Ax  (AS) 

by  partial  Gauss  Jordan  elimination,  i.e.,  by  solving  (AS)  for  pi,  pt,  ...,  pm.  Initially,  (  =  z  and 
p  =  r,  i.e.,  one  is  at  the  vertex  x  =  0.  The  vertex  specified  is  feasible  exactly  when  P  >  0  (since 
that  makes  all  the  basic  variables  non-negative  while  the  non-bask  ones  are  sero  by  choice).  In 
order  to  keep  track  of  the  value  of  the  objective  Auction,  x,  one  expresses  it  always  in  terms  of 
the  non-basic  variables, 

*  =  (A4) 

Initially,  i  =  c  and  6  =  d.  For  the  computations,  one  keepe  the  aasential  information  in  the  Ttacker 
tablean,  as  folkws: 


Pi 

-fl 

flu 

-(i 

hi 

...  i 

...  fli  ' 

Pi 

hi 

hi 

...  Pt 

Pm 

tml 

•  •  •  1«»  Pm 

X  ' 

\  71 

7l 

•••  In  f  > 

The  simplex  method  operates  on  this  tableau  in  the  following  way:  if  each  entry  in  the  final 
column  is  non-negative,  thro  one  is  currently  located  at  a  feasible  vertex.  One  wishes  to  move  to 
a  neighboring  vertex  where  the  value  of  the  objective  function  will  be  less.  If  one  of  the  entries  in 
the  final  row  is  negative,  then  one  can  increase  the  value  of  the  correspond  in  g  non-basic  variable 
and  reduce  the  sise  of  the  objective  function.  This  is  done  by  exchanging  the  non-bask  variable 
with  a  basic  variable.  This  amounts  to  solving  the  equation  representing  the  basic  variable  for  the 
non-basic  variable  and  then  substituting  the  resulting  expression  into  all  of  the  other  equations. 

Suppose  one  wishes  to  exchange  the  variable  in  the  Jfc-th  row  with  the  variable  in  the  1-th 
column.  The  folkwing  pivot  roles  can  be  derived  to  carry  out  this  process.  Let  T  represent  the 
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entire  tableau,  and  tt]  its  t,j  entry.  Then 


tu  -*  l/*w 

tu  *  tn/tici  t  =  1, ...,  k  —  1,  k  ■+■  1, m  +  1  <^6» 

tkj  -tkj/tki  ’  j  =  1, ...,/- 1,/+ l,...,n+ 1 

t%j  —*  tij  —  tutifj/tki 

In  order  to  maintain  feasibility,  one  must  choose  the  pivot  row  carefully.  This  is  done  by  selecting 
the  row  k  by  the  following  rule: 

Choose  k  so  that  _  mn  ( i  ti,  >  o\ .  (A7) 

tki  *  l  Ui  ) 

If  more  than  one  row  satisfies  this,  then  any  of  them  may  be  chosen.  One  performs  exchanges  until 
the  final  row  of  the  tableau  is  non-negative.  One  then  has  an  optimal  solution. 

Of  course,  if  the  rightmost  column  has  some  negative  entries,  then  there  is  some  extra  work 
involved  before  one  can  start  this  optimization  process;  one  must  first  find  a  feasible  vertex.  For¬ 
tunately,  the  early  pioneers  of  linear  programming  noted  that  the  simplex  method  can  even  be 
used  to  tackle  thiB  problem.  After  all,  it  is  a  procedure  which  transforms  a  tableau  into  one  whose 
last  row  is  non-negative,  while  preserving  the  non-negativity  of  the  left  column.  So,  one  simply 
considers  the  dual  simplex  method,  in  which  the  roles  of  basic  and  non-basic  variables  are  re¬ 
versed.  A  new  objective  function  (0  works  well  because  it  eliminates  the  necessity  of  adding  a  row 
to  the  tableau)  is  added  to  the  problem.  This  objective  function  should  have  the  property  that  its 
coefficients  in  terms  of  the  current  non-basic  variables  are  non-negative,  but  no  other  restrictions 
need  to  be  placed  on  its  selection.  Now  one  has  a  tableau  whose  last  row  is  non-negative,  but 
whose  last  column  does  not  share  this  property.  One  can  apply  the  dual  simplex  method  to  this 
tableau  to  get  one  in  which  both  the  last  row  and  the  last  column  are  non-negative.  The  dual 
simplex  method  proceeds  just  as  the  simplex  method  proceeds,  except  that  the  row  and  column 
pivot  rules  are  interchanged,  the  former  pivot  column  selection  rule  is  now  used  to  choose  the  pivot 
row,  and  the  column  /  is  chosen  by  the  following  rule: 

Choose  /  so  that  -  *JS±LL  =  nun  tkj  <  0)  .  (A8) 

tki  3  l  tkj  ) 

One  continues  exchanging  basic  and  non-basic  variables  until  the  rightmost  column  is  non-negative. 
At  this  point,  one  has  a  feasible  vertex,  and  the  new  objective  function  can  be  replaced  by  the 
original,  which,  of  course,  must  be  expressed  in  terms  of  the  current  non-basic  variables,  something 
best  accomplished  by  keeping  track  of  it  all  along.  One  can  now  proceed  with  the  optimization 
process  as  before. 

This  entire  procedure,  as  well  as  many  of  its  alternatives,  is  described  in  Dantzig  [1963]. 
Mangas arian  [1978]  discusses  the  Tucker  tableau  and  many  of  its  applications. 
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